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STIR  PROJECT: 

TRANSIENT  SOIL  DENSITY  IMPACTS  LAND  SURFACE  CHARACTERISTICS  AND 

CHARACTERIZATION 


Y.  Kojima,  J.L.  Heitman  (PI),  and  R.  Horton 


ABSTRACT 

Soil  density  is  commonly  treated  as  static  in  studies  on  land  surface  property  dynamics. 
Magnitudes  of  errors  associated  with  this  assumption  are  largely  unknown.  Objectives  of  this 
preliminary  investigation  were  to:  i)  quantify  effects  of  soil  density  variation  on  soil  properties, 
and  ii)  evaluate  impact  of  changing  soil  density  on  surface  energy  balance  and  heat  and  water 
transfer.  Six  soil  properties  were  evaluated  over  a  range  of  soil  densities,  using  a  combination  of 
ten  modeling  approaches.  Thermal  conductivity,  water  characteristics,  hydraulic  conductivity, 
and  vapor  diffusivity  were  most  sensitive;  these  properties  changed  by  fractions  greater  than 
associated  change  in  density  (i.e.,  10%  change  in  density  led  to  >10%  change  in  property). 
Subsequently,  three  field  seasons  were  simulated  with  a  numerical  model  (HYDRUS-1D)  for  a 
range  of  bulk  densities.  Among  the  surface  energy  balance  terms,  ground  heat  flux  and  latent 
heat  flux  were  most  sensitive  to  bulk  density.  Surface  soil  temperature  variation  increased  in 
with  low  bulk  densities  but  subsurface  temperature  variation  decreased.  Surface  water  content 
varied  with  bulk  density  but  effects  mostly  disappeared  in  the  subsurface.  Results  demonstrate 
significance  of  transient  density  on  surface  conditions  and  point  to  need  for  continued  evaluation 
of  impacts  with  improved  measurements  and  modeling. 
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INTRODUCTION 

Surface  soil  is  a  complex,  dynamic  interface  which  dictates  mass  and  energy  transfer  between 
land  and  atmosphere,  and  determines  water  flow  and  partitioning  in  the  hydrological  cycle.  Its 
properties  are  considered  dynamic  because  they  are  controlled  in  part  by  soil  water  content, 
which  can  change  quickly  with  wetting  events  or  slowly  over  sustained  periods  of  drainage,  plant 
uptake,  and  evaporative  drying.  Filling  and  emptying  of  water  in  soil  pore  space  alters  soil 
hydraulic  and  thennal  properties.  Because  understanding  soil-water  controlled  properties  is 
critical  for  modeling  and  interpreting  broader  hydrologic  and  environmental  processes, 
tremendous  effort  has  been  expended  to  develop  soil  water  sensor  technologies  and  monitoring 
networks  (Robinson  et  al.,  2008;  Ochsner  et  ah,  2013).  This  work  has  led  to  new  understanding 
of  soil  property  dynamics,  and  for  potentially  even  greater  understanding,  as  these  measurements 
are  coupled  with  remote  sensing  to  extend  measurement  footprints  (Albergel  et  al.,  2012).  Yet,  in 
all  of  these  efforts  there  remain  fundamental  questions  that  have  not  been  addressed.  An 
elementary  and  ubiquitous  assumption  in  hydrologic  studies  considering  dynamic  soil  surface 
properties  is  that  soil  density  is  static.  We  know,  in  fact,  that  this  is  not  the  case. 

Consequences  associated  with  this  assumption  are  largely  unknown,  but  are  likely  critical  (cf. 
Arya  and  Paris,  1981;  Moldrup  et  al.,  2000;  Ochsner  et  al.,  2001).  Large  areas  of  the  land  surface 
undergo  significant  changes  in  surface  soil  density  through  annual  cycles  of  disturbance 
associated  with  agriculture  (Strudley  et  al.,  2008;  Logsdon,  2012;  Liu  et  al.,  2014).  Freeze-thaw 
processes  alter  surface  density  and  arrangement  seasonally  (Staricka  and  Benoit,  1995).  Shrink- 
swell  processes,  erosion,  and  deposition  alter  surface  soil  density  and  arrangement  episodically 
(Timm  et  al.,  2006).  Unfortunately,  due  to  historical,  practical  limitations  in  our  ability  to 
continuously  quantify  soil  density-derived  effects,  this  limitation  remains  mostly  unaddressed  as 
a  dynamic  factor  in  land  surface  characterization,  and  the  magnitudes  of  any  associated  errors  are 
unknown  to  scientists  and  engineers  working  on  a  multitude  of  related  investigations. 

Our  general  goal  is  to  examine  the  impact  of  transient  (i.e.,  dynamic)  soil  density  on  land  surface 
characteristics  and  characterization.  Moving  forward  toward  this  goal  likely  requires  both 
extensive  measurement  and  modeling  efforts.  The  objectives  for  this  STIR  project  were  focused 
toward  making  initial  progress  in  work  aimed  at  this  longer-term  goal. 

Specific  objectives  for  this  project  were  to: 

1)  Quantify  effects  of  soil  density  variation  on  fundamental  soil  properties. 

2)  Evaluate  impact  of  changing  soil  density  and  associated  properties  on  surface  energy  balance 
and  coupled  heat  and  water  transfer  in  soils. 

We  first  modeled  a  series  of  gennane  soil  properties:  volumetric  heat  capacity,  thennal 
conductivity,  soil  thermal  diffusivity,  water  retention  characteristics,  hydraulic  conductivity,  and 
vapor  diffusivity  using  soil  property  models  from  the  literature  that  included  the  capacity  to 
incorporate  bulk  density  dependence.  We  then  used  a  subset  of  these  properties  in  a  numerical 
model  to  examine  expected  variation  in  surface  energy  balance  terms,  soil  water  content,  and  soil 
temperature  associated  with  bulk  density  variation  as  a  case  study.  In  the  following,  methods  and 
results  are  summarized  together  by  analysis.  Detailed  interpretation  is  currently  being  pursued 
for  a  forthcoming  manuscript. 
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PROPERTY  MODELING  (Objective  1) 

We  based  our  analysis  on  the  range  of  bulk  densities  (pb)  observed  in  several  previous  field 
studies.  Bemdt  and  Coughlam  (1976)  investigated  pb  changes  associated  with  shrink-swell  of  a 
clay  soil  and  reported  values  ranging  from  1.04  to  1.37  Mg  m  3  (32%  variation)  with  wetting¬ 
drying  cycles.  Kays  et  al.  (1985)  showed,  for  a  clay  loam  soil,  that  pb  change  associated  with 
freeze-thaw  cycles  was  from  1.28  to  1.18  Mg  m'3  (8%  variation).  Logsdon  (2012)  observed  pb 
variation,  through  periodic  sampling  over  5  yr  in  a  clay  loam,  of  1.25  to  1.40  Mg  m  (12% 
variation).  Liu  et  al.  (2014)  observed  pb  variations  due  to  settling  after  tillage  and  reported  values 
for  silt  loam  changing  from  1.00  to  1.35  Mg  m'  (35%  change),  and  values  for  a  sandy  loam 
changed  from  1.10  to  1.35  Mg  m  3  (23%  change). 

We  tested  three  soils  in  our  analysis  (Table  1):  silt  loam  and  sandy  loam  with  data  from  Liu  et  al. 
(2014),  and  clay  loam  from  Logsdon  et  al.  (2012)  using  properties  of  a  Webster  soil  series. 

Soil  Thermal  Properties 

Impact  of  pb  on  soil  volumetric  heat  capacity  C  was  evaluated  with  the  de  Vries  (1963)  model, 


C  =  csPb  +  CL0  [1] 

where  cs  is  specific  heat  of  soil  solid  (J  kg  1  °C_1),  Cl  is  volumetric  heat  capacity  of  liquid  water 
(4,184,000  J  m3  K  ')  and  9  is  volumetric  water  content  (m3  m3).  The  values  for  cs  were 
calculated  based  on  particle  size  distribution  and  soil  organic  matter  content  (SOM)  as  described 
by  de  Vries  (1963).  The  C  values  were  calculated  for  9  =  0.10  m  m  to  0.40  m  m'  and  with  pb 
=  1.10  Mg  m  3  to  1 .40  Mg  m  3. 

As  would  be  expected  from  the  model,  increasing  pb  resulted  in  an  increase  of  C,  and  the  rate  of 
increase  was  constant  over  the  range  of  0  (Fig.  1).  If  we  treat  our  minimum  pb  as  the  standard  and 
use  0  =  0.25  m  m  for  evaluation,  a  100  kg  m'  increase  in  pb  results  in  77-84  kJ  m'  K"  increase 
in  C.  In  percent,  approximately  8%  increase  in  pb  resulted  in  4.0,  4.1,  and  4.0  %  increase  in  C  for 
clay  loam,  silt  loam,  and  sandy  loam,  respectively.  The  C  values  for  observed  maximum  pb  were 
106%,  110%,  and  108%  of  the  values  for  observed  minimum  pb  for  clay  loam,  silt  loam,  and 
sandy  loam,  respectively. 


Impact  of  pb  on  soil  thermal  conductivity  k  was  evaluated  with  the  de  Vries  (1963),  Campbell 
(1985),  and  Lu  et  al.  (2014)  models,  de  Vries  (1963)  provided  an  equation  to  calculate  X  based  on 
the  soil  particle  size  distribution,  pb,  and  soil  organic  matter  content: 


X  = 


X  lokixili 

X  -=oM, 


[2] 


where  N  is  the  number  of  types  of  soil  constituents,  kj  is  a  weighting  factor,  x,  is  volume  fraction, 
and  7,  is  thennal  conductivity  of  each  soil  constituent.  Empirically  detennined  values  for  k\  were 
used  in  this  study. 

Campbell  (1985)  provided  the  following  equation 
X  =  A  +  BQ  —  (A  —  D)  exp  [—  (£0)4  ]  [3] 

where  A,  B ,  D,  and  E  are  shape  factors  associated  with  soil  properties.  Empirical  parameters  can 
be  calculated  as: 
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0.57  +  1.730n  +0.939m 

A=  l-0.749,;0.499,„  -2-89-(1-9-> 
B  =  2.80 


£>  =  0.03  +  0.79,,, 


E  =  \  +  2.6mc  [7] 

where  0q  is  volume  fraction  of  quartz,  0m  is  volume  fraction  of  other  minerals,  0SOHd  is  volume 
fraction  of  soil  solid,  and  mc  is  clay  fraction.  This  model  does  not  account  for  soil  organic  matter. 
In  this  study  we  assumed  that  0q  is  equal  to  volume  fraction  of  sand,  and  0m  is  equal  to  the 
volume  fraction  of  silt  plus  clay. 

Lu  et  al.  (2014)  provided  the  following  equation 
k  =  kdiy+exp(p-0-“)  [8] 

where  Ciry  is  thermal  conductivity  of  oven  dried  soil,  and  a  and  P  are  shape  factors.  The  thermal 
conductivity  of  oven  dried  soil  can  be  estimated  from  soil  porosity  x. 

^diy  =  — 0.56x  +  0.51  [9] 

The  shape  factors  a  and  P  can  be  detennined  based  on  soil  particle  distribution  and  pb 
a  =  0.67 /cl  +  0.24  [10] 

P  =  l-97/ga  +0.00187pb  —0.00136 /saPb  -0.95  [11] 


where  fc\  is  clay  mass  fraction,  fsa  is  sand  mass  fraction. 

For  the  de  Vries  (1963)  model,  X  increases  with  increases  of  pb  and  0  (Fig.  2).  Values  of  X  were 
highest  for  sandy  loam,  and  lowest  for  silt  loam  (sandy  loam>clay  loam>  silt  loam).  When  0  = 
0.25  m3  m  3,  a  100  kg  m  3  (8%)  increase  in  pb  results  in  12.8%,  11.9%,  and  14.4%  increases  in  X 
for  clay  loam,  silt  loam,  and  sandy  loam,  respectively.  The  X  values  for  observed  maximum  pb 
were  122%,  140%,  and  131%  of  the  values  for  observed  minimum  pb  for  clay  loam,  silt  loam, 
and  sandy  loam,  respectively. 

The  Campbell  (1985)  model  is  shown  in  Fig.  3.  Trends  were  similar  to  those  from  the  de  Vries 
(1963)  model,  but  values  of  X  were  generally  larger.  Values  of  X  were  highest  for  sandy  loam, 
and  lowest  for  silt  loam  (sandy  loam>silt  loam>  clay  loam).  With  0  =  0.25  m3  m  3,  a  100  kg  m  3 
increase  in  pb  caused  11.2%,  11.0%,  and  11.8%  increases  of  X  for  clay  loam,  silt  loam,  and  sandy 
loam,  respectively.  The  X  values  for  observed  maximum  pb  were  117%,  130%,  and  125%  of  the 
values  for  observed  minimum  pb  for  clay  loam,  silt  loam,  and  sandy  loam,  respectively.  These 
percentage  increases  are  smaller  than  those  with  the  de  Vries  model,  despite  similar  changes  in  X, 
because  of  generally  higher  X  values  for  the  Campbell  model. 

3 

The  Lu  et  al.  (2014)  model  followed  similar  trends  as  the  other  models  (Fig.  4).  With  0  =  0.25  m 
m  ,  a  100  kg  m  increase  in  pb  caused  15.4%,  16.4%,  and  12.0%  increases  in  X  for  clay  loam, 
silt  loam,  and  sandy  loam,  respectively.  The  X  values  for  observed  maximum  pb  were  124%, 
146%,  and  125%  of  the  values  for  observed  minimum  pb  for  clay  loam,  silt  loam,  and  sandy 
loam,  respectively. 
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Thermal  diffusivity  k  (=)JC)  was  calculated  as  a  function  of  pb  and  0  based  on  the  de  Vries 
(1963)  model  for  C  and  Campbell  (1985)  model  for  X  (Fig.  5).  The  effect  of  pb  is  very  small  with 
0  <  0.15  m  m  for  clay  loam  and  silt  loam,  but  generally  increases  in  pb  resulted  in  increases  in 
k.  With  0  =  0.25  nr  m  ,  a  100  kg  m  3  increase  in  pb  caused  7.0%,  6.6%,  and  7.5%  increases  in  k 
for  clay  loam,  silt  loam,  and  sandy  loam,  respectively.  The  k  values  for  observed  maximum  pb 
were  111%,  118%,  and  116%  of  values  for  observed  minimum  pb  for  clay  loam,  silt  loam,  and 
sandy  loam,  respectively. 


Soil  Hydraulic  Properties 

Water  characteristics  were  first  examined  at  different  values  of  pb  and  0  using  ROSETTA,  which 
is  a  hierarchical  pedotransfer  function  (Schaap  et  ah,  2001).  Empirical  parameters  for  the  van 
Genuchten  (1980)  water  retention  model:  0S,  0r,  a,  and  n  are  output  by  ROSETTA.  ROSETTA 
also  outputs  saturated  hydraulic  conductivity  Ks.  The  van  Genuchten  (1980)  model  is 


0  =  0, +(0.-0,) 


— I  1 — 1  /  « 


1  +  a\\i‘ 


[12] 


where  vp  (m)  is  soil  water  matric  potential,  0S  and  and  0r  are  often  referred  to  as  saturated  and 
residual  water  contents,  respectively. 

Increases  in  pb  caused  decreases  in  Qr,  0S,  and  Ks  (Table  2).  Decreases  in  0S  and  Ks  are  reasonable 
but  values  for  0r  are  expected  to  increase  because  of  reduction  of  pore  size  (Assouline,  2006a). 
Figure  6  shows  resulting  soil  water  retention  curves  as  a  function  of  pb  and  v| /.  Increases  in  pb 
shift  the  water  retention  curves  downward.  A  significant  impact  of  pb  increase  is  shown  clearly  in 
the  decrease  of  saturated  water  contents.  There  are  also  relatively  large  differences  in  matric 
potential  when  soil  is  dry.  With  vp  =  -10  m,  water  content  decreased  0.002-0.008  m  m  for  a  100 
kg  m'3  increase  in  pb.  Impact  of  altering  pb  is  more  significant  in  finer  textured  soil,  i.e.,  clay 
loam. 

The  effect  of  pb  changes  on  water  characteristics  was  also  tested  with  the  model  suggested  by 
Assouline  (2006a).  Assouline  (2006a)  described  the  water  retention  curve  as 

Se  =  l-e?q)j-[a(jv|/|  1  -kLr)f}  [13] 

where  Se  is  degree  of  saturation,  Se  =  (0-Or)/(0s-0r),  a  and  p  are  fitting  parameters,  i|//  is  matric 
potential  corresponding  to  a  very  small  water  content,  0l,  which  represents  the  limit  of  the 
domain  of  interest  of  the  water  retention  curve.  For  convenience,  0r  can  be  assumed  to  equal  dL. 

Brooks  and  Corey  (1964)  suggested  the  following  expression  of  the  water  retention  curve 


Se=(vA%)CT  ,141 

Se  =1  vp  >  v|/a 

where  vpa  is  air  entry  pressure,  and  a  is  a  pore-size  distribution  index.  Assouline  showed  that  a  is 
related  to  parameters  in  Eq.  [13] 

cr  =  0.8  ls“°  837  [15] 
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8  = 


[16] 


(a^r^rxi  +  1/fa)  +  l/|v|/J 

{(oO~2/,x[r(i  +  2/fa)  —  r2(i  + 1/ fa)]}12 

Assouline  (2006a)  presented  equations  where  water  retention  curves  at  varying  pb  can  be 
estimated  when  fitting  parameters  a,  p,  and  \|/a  are  determined  with  experimental  data  at  one  pb. 
The  parameters  in  Eqs.  [  13]-[  14]  for  the  new  water  retention  curve  with  different  pb,  ac,  pc,  and 
vf/ac,  are  described  as 


ac  =  a(pbc/pb)372 

[17] 

Pc  =  P(pbc/Pb)“ 

[18] 

Vac  =  Ya(pbc/Pb)3'82 

[19] 

where  pb  and  pbc  are  original  and  new  bulk  densities,  and  co  is  defined  as 
co  =  2.3  - 1 .9  (SC/CCY0-5  [20] 

where  SC  and  CC  are  mass  fraction  of  silt  and  clay.  The  values  for  9S  and  9r  also  change  with 
changes  in  pb.  Assouline  (2006a)  presented 

=0,[(p.  -PbcMp,  -Pb)]  pi] 

e,c  =e,(Pbc/Pb)  P2] 

where  9SC  are  9rc  are  saturated  and  residual  water  content  with  the  new  pb  value,  and  ps  is  soil 
solid  density  (-2650  kg  m  3). 

In  this  study,  we  first  obtained  water  retention  parameters  for  the  van  Genuchten  (1980)  model 
from  ROSETTA  for  clay  loam,  silt  loam,  and  sandy  loam  at  pb  =  1.25  Mg  m  .  Parameters  for 
Eqs.  [  13]-[  14]  were  determined  by  fitting  data  from  these  water  retention  curves.  Based  on  these 
parameters  and  new  pb,  water  retention  curves  were  estimated  with  Eqs.  [15]-[21], 

Table  3  shows  a  subset  of  the  parameters  required  for  Eqs.  [  12]-[  14]  at  different  pb.  Estimates  of 
9r  were  plausible  by  this  approach  in  that  9r  decreases  as  pb  decreases.  Thus,  the  water  retention 
curves  with  different  pb  values  cross  one  another  (Fig.  7).  At  the  same  \p,  water  content 
sometimes  increased  and  sometimes  decreased,  depending  on  vp  and  soil  type.  For  example, 
when  vp  =  -1  m  and  with  a  199  kg  m  3  increase  in  pb,  water  content  increased  9.92  in  clay  loam, 
water  content  decreased  9.93  in  silt  loam,  and  water  content  increased  9.91-9.92  in  sandy  loam. 
When  \p  =  -10  m  and  with  a  199  kg  m  3  increase  in  pb,  water  content  increased  9.994  in  clay 
loam,  water  content  decreased  9.994  in  silt  loam,  and  water  content  decreased  9.91  in  sandy 
loam. 

Hydraulic  conductivity  K  can  be  expressed  by  the  van  Genuchten  model  (1989)  and  parameters 
provided  by  ROSETTA  as 

/f(S„)=  A'”)"]'  [23] 

where  Se  for  Eq.  [23]  is 
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1 


1—1  In 


[24] 


S.  = 


1  +  |a\|/|" 


Hydraulic  conductivity  as  a  function  of  vj/  and  pb  is  shown  in  Fig  8  (and  Ks  is  given  in  Table  2). 
As  would  be  expected,  K  consistently  decreased  with  increasing  pb. 

Assouline  (2006b)  presented  an  alternate  model  to  describe  Ks  as  a  function  of  pb,  the  value  Ksc  is 


e..o  -  e,„ ) 

2.5 

P  \ 

Va 

2 

CTcCl +cr) 

V  0s  -  0r  y 

cr(]  +  crc  ) 

[25] 


where  Ks,  0s,  0r,  vpa,  and  a  are  saturated  hydraulic  conductivity,  saturated  water  content,  residual 
water  content,  air  entry  pressure,  and  pore  size  distribution  index  at  standard  pb,  and  0SC,  0rc,  \|/ac, 
and  <7C  are  saturated  water  content,  residual  water  content,  air  entry  pressure,  and  pore  size 
distribution  index  at  the  new  pb.  Unsaturated  hydraulic  conductivity  K(SC)  for  a  variety  of  pb 
were  estimated  with  the  Mualem  (1976)  and  Brooks  and  Corey  (1964)  model 

K(Se  )  =  KsSei2+2"5cT)/cT  [26] 

Values  of  Ks  for  different  pb  are  shown  in  Table  3.  Note  that  Ks  at  1.25  Mg  m  was  used  as  the 
standard  for  modification  at  other  pb  according  to  Eq.  [26].  Figure  9  shows  K  estimated  with  Eqs. 
[25]  and  [26]  as  a  function  of  pb  and  v| /  instead  of  Se.  Relative  effects  —  K  decreases  with 
increasing  pb  —  are  similar  to  those  obtained  from  ROSETTA,  despite  differences  in  the  water 
characteristics  discussed  earlier  because  the  Mualem  model  treats  residual  water  content  as 
immobile. 

Vapor  diffusivity  Dv  in  soil  can  be  described  as  (Saito  et  ah,  2006) 

Dv=zQaDa  [27] 


where  x  is  a  tortuosity  factor,  0a  is  air  filled  porosity,  and  Da  is  water  vapor  diffusivity  in  air.  The 
tortuosity  factor  can  be  described  as  (Millington  and  Quirk,  1961) 


T 


M  7/3 


[28] 


Since  0a  and  0S  are  simple  functions  of  pb  and  soil  water  content,  there  is  no  influence  of  different 
soil  type  on  Dv.  Values  for  Dv  decrease  with  increasing  pb  because  of  the  associated  decrease  in 
0a  (Fig.  10).  When  0  is  0.25  nr  m'3,  a  100  kg  m  3  increase  of  pb  caused  28.6%,  23.3%,  and  24.9% 
decrease  in  Dv  for  clay  loam,  silt  loam,  and  sandy  loam,  respectively.  The  Dv  values  for  observed 
maximum  pb  were  58.8%,  47.2%,  and  53.6%  of  the  values  for  observed  minimum  pb  for  clay 
loam,  silt  loam,  and  sandy  loam,  respectively. 

Key  Findings  from  Property  Modeling 

Six  soil  thermal  and  hydraulic  properties  were  evaluated  for  three  soil  textures  over  a  realistic 
range  in  transient  field  soil  bulk  density,  using  a  combination  of  ten  models/modeling  approaches 
available  from  the  literature.  The  properties  that  appeared  to  be  most  sensitive  to  bulk  density  are 
as  follows: 


•  Thermal  conductivity  -  change  of  <10%  in  bulk  density  led  to  11-16%  change  in  thermal 
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conductivity. 

•  Water  characteristics  -  25%  change  in  bulk  density  led  to  20-25%  change  in  residual  and 
saturated  water  contents,  with  changes  occurring  in  opposite  directions  (i.e.,  larger  residual 
water  content  and  smaller  saturated  water  content). 

•  Saturated  hydraulic  conductivity  -  values  for  saturated  hydraulic  conductivity  typically 
change  by  an  order  of  magnitude  over  the  range  of  transient  field  bulk  density. 

•  Vapor  diffusivity  -  change  of  <10%  in  bulk  density  led  to  23-29%  change  in  diffusivity. 


NUMERICAL  SIMULATIONS  (Objective  2) 

Simulations  were  performed  with  the  HYDRUS-1D  software  package  (Simunek  et  ah,  2009)  to 
evaluate  impacts  of  change  in  bulk  density  on  surface  energy  balance  and  soil  heat  and  water 
transfer.  Four  soil  profiles  (A,  B,  C,  and  D)  were  used  in  the  simulations,  each  approximately 
representing  a  soil  with  silt  loam  texture.  The  soil  profiles  have  two  layers  (Fig.  11),  one 
represents  a  disturbed  soil  layer  (0-0.225  m  depth)  which  has  pb  =  1.3  (A),  1.2  (B),  1.4  (C)  or  1.5 
(D)  Mg  m  3  bulk  density,  and  the  other  is  an  undisturbed  deep  soil  layer  (0.225-5  m  depth)  which 
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has  pb  =  1.4  Mg  m  .  (Thus,  profile  (C)  has  uniform  pb  throughout  the  profile.)  Node  spacing  was 
0.01  m  from  surface  through  50  cm  depth,  and  node  spacing  was  gradually  increased  to  a 
maximum  of  0.05  m  below  50  cm  depth.  Hydraulic  properties  were  expressed  with  the  Brooks 
and  Corey  (1964)  model,  and  parameters  were  obtained  with  the  Assouline  (2006a)  approach 
described  above.  Thennal  properties  were  calculated  with  the  Campbell  (1985)  model. 

Weather  data  from  an  experimental  field  near  Ames,  IA  in  2012,  2013,  and  2014  were  used  to 
detennine  surface  boundary  conditions.  Calculations  were  made  with  data  during  May-October 
in  each  year.  These  three  years  provide  differing  amounts  of  precipitation  during  May-October. 
Accumulated  precipitation  in  May-October  was  337,  524,  and  801  mm  in  2012,  2013,  and  2014, 
respectively,  i.e.,  2012  was  a  dry  year,  2013  was  intermediate,  and  2014  was  wet.  Accumulated 
solar  radiation  during  May-October  was  3781,  3422,  and  3216  MJ  m'2  in  2012,  2013,  and  2014, 
respectively.  The  dry  year  (2012)  had  greater  accumulated  solar  energy. 

The  soil  surface  boundary  condition  was  determined  by  the  calculated  surface  energy  balance 
and  the  observed  precipitation.  The  calculation  processes  are  described  in  Simunek  et  al.  (2009). 
The  bottom  boundary  conditions  were  free  drainage  for  water  transport  and  zero  gradient  for  heat 
transport.  The  initial  condition  for  water  content  was  0.25  nr  m  3  for  all  depths  and  the  initial 
condition  for  temperature  was  20°C  at  all  depths. 

Surface  Energy  Balance 

Across  all  three  years,  Net  radiation  was  smaller  during  the  daytime  and  larger  at  night  when  pb 
was  low  (not  shown).  The  relatively  small  differences  were  likely  associated  with  differences  in 
surface  albedo  and  longwave  radiation  from  soil  surface.  Ground  heat  flux  showed  relatively 
large  differences  with  pb  variation,  particularly  for  dynamic  fluxes  within  a  given  day  (Fig.  12). 
Accumulated  differences  on  an  annual  basis  were  relatively  small  (Table  4).  In  2013  and  2014, 
ground  heat  flux  was  relatively  small  at  pb  =  1 .2  Mg  m  3,  and  generally  increased  with  pb  (Table 
4).  This  may  be  associated  with  greater  thermal  conductivity  with  larger  pb.  However,  trends 
differed  in  2012  when  conditions  were  driest.  In  most  cases  smaller  pb  produced  larger  latent 
heat  flux  (not  shown).  Accumulated  latent  heat  flux  (calculated  as  evaporation  depth)  was  the 
highest  with  the  lowest  pb  in  each  simulated  year  (Table  4).  This  trend  likely  corresponds  with 
increased  storage  of  water  available  for  evaporation  from  changes  to  the  water  characteristics 
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and  with  greater  vapor  diffusivity  at  low  pb-  Based  on  these  differences,  surface  energy 
partitioning  shifted  toward  a  relatively  greater  proportion  of  available  energy  partitioning  to 
sensible  heat  flux  when  pb  was  largest. 

Soil  Heat  and  Water  Dynamics 

Soil  temperature  at  the  5  cm  soil  depth  generally  showed  greatest  daily  variation  with  low  pb 
(Fig.  13).  Differences  in  temperature  at  maximum  and  minimum  were  typically  on  the  order  of  1 
°C  with  low  pb  having  both  the  largest  maximum  and  smallest  minimum  (i.e.,  difference  in 
variation  of  2  °C).  At  the  30  cm  soil  depth,  where  pb  was  the  same  for  each  simulated  profile, 
surface  pb  also  influenced  observed  temperatures  (Fig.  14).  However,  in  this  case  the  trend  was 
opposite  that  observed  at  the  surface.  At  the  30  cm  depth,  daily  temperature  variation  increased 
with  high  surface  pb.  In  this  case,  the  surface  layer  with  low  pb,  and  thus  low  thermal 
conductivity,  acts  as  insulation,  muting  temperature  variation  in  the  subsurface.  On  a  seasonal 
basis,  high  surface  pb  results  in  earlier  warming  in  the  summer  and  earlier  cooling  in  the  fall  (Fig. 
15). 

Soil  water  content  at  the  5  cm  depth  was  generally  drier  at  low  pb  (Fig.  16).  This  result  is  likely  a 
combination  of  both  more  rapid  drainage  during  rainfall  events,  and  lower  residual  water  content 
retained.  During  a  typical  drying  event,  simulated  water  content  at  low  pb  was  about  0.02  nr  m  3 
lower  than  at  the  largest  pb  (Fig.  17).  At  the  30  cm  depth,  differences  between  profiles  with 
different  pb  were  generally  small  (not  shown). 

Key  Findings  from  Numerical  Simulation 

Three  seasons  with  differing  surface  conditions  (rainfall,  solar  radiation)  were  simulated  with  a 
numerical  model  for  a  range  of  bulk  density  conditions.  Main  findings  were  that  as  bulk  density 
increased: 

•  Ground  heat  flux  increased  by  as  much  as  25%  on  an  annual  basis,  though  effects  varied  by 
year. 

•  Evaporation  rate  (latent  heat  flux)  decreased  by  as  much  as  7-8%  on  an  annual  basis. 

•  Surface  layer  temperature  variation  decreased  -  differences  in  variation  at  the  5  cm  depth 
were  on  the  order  of  2  °C. 

•  Subsurface  layer  temperature  variation  increased  -  even  at  30  cm  depth,  the  effect  was  on  the 
order  of  1  °C. 

•  Surface  soil  water  content  increased  by  about  0.02  m  m  during  typical  drying  events. 
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Table  1.  Soil  particle  size  distribution,  soil  organic  matter  content  (SOM),  and 
observed  minimum  and  maximum  values  for  soil  bulk  density. _ 


Particle  size  distribution 

SOM 

Bulk  density 

Texture 

Sand 

Silt 

Clay 

Min 

Max 

kg  kg  1 

kg  m 3 

Clay  loam 

21 

47 

32 

0.07 

1250 

1400 

Silt  loam 

17 

62 

21 

0.01 

1100 

1350 

Sandy  loam 

53 

38 

9 

0.01 

1150 

1350 

Table  2.  Empirical  parameters  for  the  van  Genuchten  (1980)  model 
output  with  ROSETTA  (Schaap  et  al.,  2001)  as  a  function  of  soil 
bulk  density. 

Parameter 

Bulk  density 

Mg  m 3 

1.1 

1.25 

1.4 

Clay  loam 

0,- 

0.091 

0.088 

0.084 

0, 

0.522 

0.480 

0.441 

Ks  (cm  d'1) 

53.2 

22.5 

9.6 

Silt  loam 

0,- 

0.078 

0.074 

0.071 

0, 

0.495 

0.455 

0.419 

Ks  (cm  d'1) 

70.5 

33.7 

15.8 

Sandy  loam 

0,- 

0.045 

0.042 

0.040 

0, 

0.442 

0.432 

0.377 

Ks  (cm  d'1) 

121 

68.8 

40.1 

11 


Table  3.  Empirical  parameters  for  the  Assouline  (2006)  approach  as 
a  function  of  soil  bulk  density.  Parameters  at  bulk  density  =  1.25 
Mg  m-3  were  estimated  with  ROSETTA  (Schaap  et  al.,  2001). 


Parameter 

Bulk  density 

Mg  m 3 

1.1 

1.25 

1.4 

Clay  loam 

0,- 

0.080 

0.088 

0.102 

0, 

0.531 

0.480 

0.428 

(cm) 

52.2 

85.0 

131 

a 

0.47 

0.52 

0.57 

Ks  (cm  d'1) 

76.1 

22.5 

Silt  loam 

6.82 

0,- 

0.067 

0.074 

0.085 

9.s 

0.504 

0.455 

0.406 

(cm) 

90.3 

147 

227 

a 

0.546 

0.63 

0.72 

Ks  (cm  d'1) 

106 

33.7 

10.8 

Sandy  loam 

0,- 

0.040 

0.042 

0.051 

0, 

0.452 

0.408 

0.365 

H 'a  (cm) 

40.4 

65.8 

102 

a 

0.42 

0.51 

0.59 

Ks  (cm  d'1) 

193 

68.8 

24.6 

Table  4.  Accumulated  ground  heat  flux  (positive  downward)  and 
evaporation  as  a  function  of  soil  bulk  density  for  May  to  October  in 

each  simulation  year. 

Soil  bulk  density  (surface) 

Mg  m  3 

1.1 

1.3 

1.4 

1.5 

Year 

Cumulative  ground  heat  flux 

MJ  m 3 

2012 

13.4 

13.3 

13.1 

13.2 

2013 

6.0 

5.7 

6.1 

6.4 

2014 

-4.6 

-4.1 

-3.7 

-3.6 

Cumulative  evaporation 

mm 

2012 

390 

380 

371 

363 

2013 

576 

569 

563 

557 

2014 

858 

854 

852 

849 
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Figure  1.  Volumetric  heat  capacity  as  a  function  of  bulk  density  and  volumetric  water  content. 
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Figure  2.  Thermal  conductivity  as  a  function  of  bulk  density  and  volumetric  water  content  with  de 
Vies  (1963)  model. 
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Figure.  3.  Thermal  conductivity  as  a  function  of  bulk  density  and  volumetric  water  content  with 

Campbell  (1985)  model. 
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Figure  4.  Thermal  conductivity  as  a  function  of  bulk  density  and  volumetric  water  content  with  Lu 

et  al.  (2014)  model. 
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Figure  5.  Thermal  diffusivity  as  a  function  of  bulk  density  and  volumetric  water  content  based  on 
thermal  conductivity  determined  with  Campbell  (1985)  model  and  heat  capacity  determined  with 
the  de  Vries  (1963)  model. 
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Figure  6.  Water  characteristics  estimated  with  ROSETTA  (Schaap  et  al.,  2001)  for  different  values 
of  soil  bulk  density. 
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(a)  Clay  loam 


Matric  potential  (-m  H?0) 


1400 

1300 

1200 

6  1100 

10  Bulk  density  (Mg  m 


(b)  Silt  loam 


Matric  potential  (-m  H„0) 


1300 


1200 


1400 


106  Bulk  density  (Mg  rrf3) 


(c)  Sandy  loam 


Matric  potential  (-m  H20) 

Figure  7.  Water  characteristics  estimated  with  Assouline  (2006a)  and  Brooks  and  Corey  (1964) 
models  at  different  values  for  bulk  density. 
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Figure  8.  Hydraulic  conductivity  estimated  with  ROSETTA  for  different  values  of  bulk  density. 
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Figure  9.  Hydraulic  conductivity  with  Mualem  (1976)  and  Brooks  and  Corey  (1964)  models  with 
different  values  for  bulk  density. 
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Figure  10.  Vapor  diffusivity  in  soil  with  different  values  for  bulk  density. 
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Figure  11.  Soil  profiled  used  in  the  HYDRUS-1D  simulations.  A  has  a  uniform  soil  properties,  and  B 
and  C  have  disturbed  layer  and  undisturbed  layer. 
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Figure  12.  Simulated  ground  heat  flux;  positive  values  indicated  downward  heat  flux.  Top:  May 
through  October,  bottom:  day  of  year  200  to  202.  Legends  indicate  bulk  density  (rhob)  in  Mg  m'\ 
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Figure  13.  Simulated  daily  soil  temperature  at  5  cm  soil  depth.  Legend  indicates  bulk  density  (pb)  in 
Mg  m'3. 
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Figure  14.  Simulated  daily  soil  temperature  at  30  cm  soil  depth.  Legend  indicates  bulk  density  (pb) 
in  Mg  m'\ 
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Figure  15.  Simulated  seasonal  soil  temperature  at  30  cm  soil  depth.  Top:  summer,  bottom:  fall. 
Legend  indicates  bulk  density  (pb)  in  Mg  m'3. 
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Figure  16.  Simulated  seasonal  soil  water  content  at  5  cm  soil  depth.  Legend  indicates  bulk  density 
(pb)  in  Mg  m 3. 
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Figure  17.  Simulated  soil  water  content  at  5  cm  soil  depth  during  a  drying  event.  Legend  indicates 
bulk  density  (pb)  in  Mg  m'\ 
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